rm(list = ls())

library(data.table)
library(estimatr)

add_backticks = function(x) {
  paste0("`", x, "`")
}

x_lm_formula = function(x) {
  paste(add_backticks(x), collapse = " + ")
}

load('./data/panel_month_dummies.RData')

panel[,losing_income_25:= ifelse(workers_losing_income_share > quantile(workers_losing_income_share, 0.25),1,0)]
panel[,losing_income_75:= ifelse(workers_losing_income_share > quantile(workers_losing_income_share, 0.75),1,0)]
panel[,covid_econ_25:=covid2*losing_income_25]
panel[,covid_econ_75:=covid2*losing_income_75]

month_cols <-  grep('month_2', colnames(panel), value=TRUE)

panel[,less_college:=(100-educ_diplBA)]
panel[,(paste('educ',month_cols,sep = '_')):= lapply(.SD, function(x) 
  x * panel[['less_college']] ), .SDcols = month_cols]

panel[,(paste('foreign',month_cols,sep = '_')):= lapply(.SD, function(x) 
  x * panel[['foreign_pop']] ), .SDcols = month_cols]


educ_cols <-  grep('educ_month', colnames(panel), value=TRUE)
educ_cols <- educ_cols[-157]
foreign_cols <- grep('foreign_month', colnames(panel), value=TRUE)
foreign_cols <- foreign_cols[-157]

panel[,time:= as.numeric(round((month-as.Date('2020-01-01',format='%Y-%m-%d'))/(365.25/12)))]

province_cols <- grep('code_province_', colnames(panel), value=TRUE)
panel[,(paste('trend',province_cols,sep = '_')):= lapply(.SD, function(x) 
  x * panel[['time']] ), .SDcols = province_cols]
trend_cols <- grep('trend_code_province', colnames(panel), value=TRUE)
trend_cols <- trend_cols[-1]

panel[,trend_econ:=time*losing_income_d]
panel[,trend_econ_25:=time*losing_income_25]
panel[,trend_econ_75:=time*losing_income_75]

flex_controls <- c(educ_cols, foreign_cols, trend_cols)

formula_flex_controls_econ <- as.formula(paste("hc_pc_asians ~ covid_econ + trend_econ +", x_lm_formula(flex_controls)))
formula_flex_controls_econ_25 <- as.formula(paste("hc_pc_asians ~ covid_econ_25 + trend_econ_25 +", x_lm_formula(flex_controls)))
formula_flex_controls_econ_75 <- as.formula(paste("hc_pc_asians ~ covid_econ_75 + trend_econ_75 +", x_lm_formula(flex_controls)))

out_econ_con_trend_all <- lm_robust(formula = formula_flex_controls_econ,
                                    data=panel, se_type='stata',
                                    clusters = panelvar, fixed_effects = ~ panelvar + month)

out_econ_25_con_trend_all <- lm_robust(formula = formula_flex_controls_econ_25,
                                       data=panel, se_type='stata',
                                       clusters = panelvar, fixed_effects = ~ panelvar + month)

out_econ_50_con_trend_all <- lm_robust(formula = formula_flex_controls_econ_75,
                                       data=panel, se_type='stata',
                                       clusters = panelvar, fixed_effects = ~ panelvar + month)


save(out_econ_con_trend_all, out_econ_25_con_trend_all, out_econ_50_con_trend_all, file = './output/out_econ_quantiles.RData')